/*
    <one line to give the library's name and an idea of what it does.>
    Copyright (C) 2013  Guillaume <email>

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/


#ifndef MSBAYSIANSCALARMODEL_H
#define MSBAYSIANSCALARMODEL_H

#include <msScalarFunction.h>
#include <msDiracs.h>
#include <uqSamplerPdf.h>
#include "msUtilitiesBayes.h"

#include <uqGslMatrix.h>
#include <uqStatisticalInverseProblem.h>
#include <uqStatisticalForwardProblem.h>
#include <uqInfoTheory.h>


namespace atomism
{
double likelihoodRoutine_(const uqGslVectorClass& paramValues,
			  const uqGslVectorClass* paramDirection, 
			  const void*  Data,
			  uqGslVectorClass*  gradVector,
			  uqGslMatrixClass*  hessianMatrix,
			  uqGslVectorClass* hessianEffect);

class uqBayesianScalarModel : public msPhysicalInterface
{ 
 friend double likelihoodRoutine_(const uqGslVectorClass& paramValues,
				  const uqGslVectorClass* paramDirection, 
				  const void*  Data,
				  uqGslVectorClass*  gradVector,
				  uqGslMatrixClass*  hessianMatrix,
				  uqGslVectorClass* hessianEffect);
 
 class uqBayesianScalarModelParams : public msParamsManager
       {
          static msRegistrar Registrar;
          static msRegister* createInstance() { return new uqBayesianScalarModelParams; }
	 
	  public:
      void initialize() { msParamsManager::initialize(); 	    
			  addParameter( "ChainLength" ,     "1000" , "Default chain length" ,    msUnit() );
			  addParameter( "LastChainLength" , "1000" , "Last level chain length" , msUnit() );
			  addParameter( "DisplayStepRatio" , "0.1" , "Display inference result ratio w/ChainLength" , msUnit() );
			  } 

      static boost::shared_ptr<uqBayesianScalarModelParams> New()
          { boost::shared_ptr<uqBayesianScalarModelParams> T( new uqBayesianScalarModelParams );
	    T->initialize(); T->update(); 
	    return T;
	    }

      uqBayesianScalarModelParams() : msParamsManager() {};
	}; 
protected:	
   
      static msRegistrar Registrar;
      static bool  isBayesianScalarModelRegisteredInPython;
      static msRegister* createInstance() { return new uqBayesianScalarModel; }
      
      void registryInPython();
      
      uqBayesianScalarModel(); 
      
public:
  
     void update()  { msTreeMapper::update();
                      };
		     
     void initialize() {  msPhysicalInterface::initialize();
                          declareChild<msGeneralizedCoordinates>(UqParameters,
								 msGeneralizedCoordinates::New(getUnits()),
								 "UqParameters");
			  declareChild<msDiracs>(Database,
						 msDiracs::New(""),
						 "Database");
			  declareChild<msScalarFunction>(Model,
						         msScalarFunction::New(""),
						         "Model");			  
			  declareChild<uqSamplerPdf>(PosteriorPDF,
						     uqSamplerPdf::New(),
						     "PosteriorPDF");
 			  PosteriorPDF->setParameters(UqParameters.getSharedPtr());
			  };
			    
     static boost::shared_ptr<uqBayesianScalarModel> New(){
             
               boost::shared_ptr<uqBayesianScalarModel> T( new uqBayesianScalarModel );
	       T->initialize();
	       T->setUnits(msUnitsManager::New());
	       T->setParameters( uqBayesianScalarModelParams::New() );
	       T->update();
	       return T;
	       }
     
     virtual boost::shared_ptr<msTreeMapper> clone() { 
            boost::shared_ptr<uqBayesianScalarModel> clone = New(); 
	    clone->msTreeMapper::operator=(*this);					      
	    return clone;
            }    
     
     boost::shared_ptr<msTreeMapper> setDatabase(boost::shared_ptr<msDiracs> database);
     
     boost::shared_ptr<msTreeMapper> setParameter(string name,double min,double max);
     
     boost::shared_ptr<msTreeMapper> setModel(boost::shared_ptr<msScalarFunction> model);
     
     boost::shared_ptr<msTreeMapper> setSampleFromJointDistribution(size_t i);
      
     boost::shared_ptr<msTreeMapper> updateModelParameters();
         

     boost::shared_ptr<msGeneralizedCoordinates> getUqParamaters() 	{ return UqParameters.getSharedPtr(); };

     boost::shared_ptr<uqSamplerPdf> getPosteriorPDF() 			{ return PosteriorPDF.getSharedPtr(); };

     void infer();
     
     ostream& print(ostream& out) const;
     
   
private:
     
     msChild<msGeneralizedCoordinates> 	UqParameters; 
     msChild<msScalarFunction> 		Model;
     msChild<msDiracs> 			Database;
     msChild<uqSamplerPdf> 		PosteriorPDF;
     

     boost::shared_ptr<uqFullEnvironmentClass> Environment;
                     
     double LikelihoodRoutine(const uqGslVectorClass& paramValues,
			      const uqGslVectorClass* paramDirection,
                              const void* functionDataPtr,
			      uqGslVectorClass* gradVector,
			      uqGslMatrixClass* hessianMatrix, 
			      uqGslVectorClass* hessianEffect);
     
     void writeInputQueso();


     double 	ErrorMin;
     double 	LikelihoodMax;
     int 	DisplayStep;
     int 	CurrentIteration;
};

}
#endif // MSDIRACS_H
